function[SD] = chavespec(Sw,dval),
% Chave's method of making the estimate robust 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End up with a bunch of modified periodograms per segment and per 
% section in the cell array. Now Chave tells us how to average these.
% First start with the median as an initial robust LOCATION estimate
% The non-robust is just the unweighted average (the mean)
Snon=mean(Sw,2);
Sloc=median(Sw,2);
SD=zeros(size(Sloc));
nwint=size(Sw,2);
%nfft = size(Sw,2);
% I'm doing the iteration for all frequencies at the same time
% Could rewrite this for every frequency at once
% Convergence criterion for all frequencies at once based on
% relative difference wrt previous iteration in percent
relperc=1;
iter=0;
while any(abs(SD-Sloc)./Sloc*100>relperc)
  disp(sprintf('Iteration %3.3i mandated by %3.3i / %3.3i frequencies',...
	       iter+1,sum(abs(SD-Sloc)./Sloc*100>relperc),length(SD)))
  if iter>0; Sloc=SD; end
  iter=iter+1;  
  % Now come up with a "practical" but "lower efficiency" SCALE estimate
  % Calculate residuals from the location estimate
  Sres=Sw-repmat(Sloc,1,size(Sw,2));
  if strcmp(upper(dval),'MAD')
    % Scale estimate is median absolute deviation of residual from
    % the median over value expected for a chi-squared distribution
    % with 2 degrees of freedom (Chave Eq. 20 and Eq. 30)
    % This is because power spectra are the sums of squares of almost
    % normally distributed variates, and hence distributed as
    % chi-squared. The number of degrees of freedom is 2 at each
    % frequency, except for the DC and Nyquist components, which have
    % only 1 degree of freedom. The chi-squared 2-distribution is
    % equivalent to the exponential distribution and thus easy to
    % calculate. We compare the scale estimates to the chi-squared
    % 2-distribution and use this as our estimate of scale.
    Sscale=median(abs(Sres),2)/(2*asinh(1/2));
  elseif strcmp(upper(dval),'IQ')
    % Scale estimate is interquartile range of residuals over expected
    % value for chi-squared distribution with 2 degrees of freedom
    % (Chave Eq. 21 and Eq. 30)
    for index=1:nfft
      Sscale(index,1)=iqr(Sres(index,:))/(2*log(3));
    end
  end

  % Now iterate with Huber weights (Chave Eq. 26)
  % How far out are the residuals in multiples of the scale estimate?
  % Use these to construct weights - far out values are downweighted
  % This number is given by Huber and assures the efficiency of the estimate
  k=1.5;
  Wght=Sres./repmat(Sscale,1,nwint);
  buv=Wght>k;
  blo=Wght<=k;
  Wght(buv)=sqrt(k*sign(Wght(buv))./Wght(buv));
  Wght(blo)=1;

  % Construct robust estimate as weighted average of segments - in one go
  % This is the new LOCATION estimate
  SD=sum(Wght.*Sw,2)./sum(Wght,2);
  %  plot(abs(SD-Sloc)./Sloc*100); pause
end
